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ON  THE  DYNAMICS  OF  NEAR-EARTH  FLIGHT 


by  Bernard  G.  Grunebaum*  Ph.D. 

Chief  of  Technical  Staff 
Land-Air,  Inc. 

The  object  of  this  paper  is  to  develop  equations  for  the  near-earth 
flight  of  missiles  and  satellites,  vvhich  permit  to  obtain  accurate  numer¬ 
ical  results.  Any  differential  equation  that  defied  analytic  solution 
had  to  be  of  the  quasi-1 inear  first  order  type  so  the  numerical  integration 
could  be  performed  with  close  error  bounds.  The  use  of  experimental  para¬ 
meters  was  avoided  as  much  as  possible  since,  at  the  present  state  of  the 
art,  the  measurement  of  such  parameters  is  rather  inaccurate  in  most  cases 
and  frequently  the  error  which  such  parameters  pretend  to  correct  is 
strongly  suspected  to  be  smaller  than  the  error  derived  from  the  measure¬ 
ments.  Both  free  flight  and  flight  under  thrust  had  to  be  considered  and 
the  equations  were  required  to  be  such  as  to  enable  us  to  analyse  numer¬ 
ically  and  separately  the  effects  of  different  physical  causes,  such  as 
earth  oblateness,  air  drag,  and  effect  of  earth  rotation  on  the  latter. 

All  analytical  data  were  assumed  to  be  determinable  from  Doppler  Radar 
and/or  Transit  equipment  and  the  algorithms  had  to  be  such  as  to  achieve 
a  reasonable  compromise  between  accuracy  and  speed. 

These  objectives  have  been  met,  at  least  from  an  analytical  stand¬ 
point.  The  paper  is  nevertheless  incomplete,  since  numerical  results  still 
have  not  been  obtained.  The  conclusion  of  this  report  list  certain  items 
of  programing  necessary  to  complete  the  study.  The  work  has  begun,  and 
when  numerical  results  have  been  obtained,  a  concluding  report  will  be 
made. 

The  writer  wishes  to  acknowledge  valuable  help  received  from  Mr.  M.  R. 
Claasen  in  verifying  some  of  the  algebraic  developments  and  from  Mr.  R. 
Bergman  in  the  form  of  suggestions  regarding  apsidal  advance. 
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Basic  Principles  of  Mechanics 

For  reasons  of  clearness,  we  shall  recall  some  of  the  basic  equations  of 
classical  dynamics.  The  following  definitions  will  be  used: 


x,y,z  =  orthogonal,  right  handed,  geocentric  coordinate  system,  x,y  is  assumed 
to  be  the  equatorial  plane  and  z  the  polar  axis  pointing  North.  The 
system  is  assumed  inertial,  i.e.  the  earth  rotates  around  the  z  axis 
in  the  x,y  direction, 
r  =  vector,  defining  a  point  P. 

P  “  projection  of  r  on  x,y  plane. 

«p  =  angle  p,  r,  positive  in  the  p,  z  direction, 

ft  =  angle  x,  p,  positive  in  the  x,y  direction. 

NBT  =  orthogonal,  right  handed,  unit  vector  system;  N  is  the  unit  vector  of 
r,  T  is  perpendicular  to  r  in  the  r,  z  plane  in  the  direction  of  in¬ 
creasing  9,  6  is  perpendicular  to  NT  in  the  direction  of  increasing  G. 
The  plane  BT  is  called  the  local  horizontal  by  definition. 

V  =  velocity  vector  of  P. 

a  =  projection  of  v  on  local  horizontal 

a  =  angle  between  v  and  local  horizontal  plane;  positive  if  v  is  on  the 

same  side  as  N.  _  . 

p  =  angle  between  a  and  B  measured  in  the  B,T  direction;  i.e.  if  *  T, 
then  P  =  j. 


Primes  indicate  differentiation  with  respect  to  time.  We  have  the 
following  basic  relations: 


V 

sin 

a  =  r 

/ 

(1) 

V 

cos 

a  cos 

P  = 

r9' 

cos 

V 

cos 

a  sin 

II 

r©' 

(2) 

a 

=  V 

cos  a 

N' 

'  = 

(p'T  +  cos  ' 

ip  9' 

'B 

(3) 

T' 

r  _ 

-f'N  - 

sin 

<p  9' 

'B 

B' 

'  3 

'  sin 

9  T 

-  9' 

'  cos  f  N 
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(4)  V*  =  r '*  +  cos*  <p. 


Furthermore,  a  can  be  easily  verified 

(5)  r*  -  r  +  rq)*T  +  r0*  cos  <p  B 
and 

(6)  r*=  (r*-  -  r0  cos*  <p)N  +  (2r*fp*  +  rq)' +  r0^*  sin  q)  cos  q))T 

+  (2r  ^0^  cos  q)  +  r0‘'  cos  q>  -  2r0  sin  q>)B 

Proofs  for  the  preceding  expressions  are  elementary  and  will  be  omitted. 

i,n^  vtrUUas 

studying  the  movement  of  a  point  in  a  field  of  forces  usually  yields 
non«l inear  differential  equations  which,  except  in  simple  and  rather  un> 
realistic  cases,  defy  analytical  integration.  Moreover,  they  are  usually 
not  of  the  quasi>l inear  first  order  type,  and  their  numerical  integration 
can  be  lengthy  and/or  of  unknown  accuracy.  It  is  well  known  that  this  is 
generally  the  case  when  the  variables  and  coordinates  of  the  previous 

paragraph  are  used  to  study  the  movement  of  a  missile  or  a  satellite. 

Therefore,  our  first  problem  is  to  choose  a  set  of  variables  and  a  system 
of  coordinates  that  will  yield  suitable  quasi-1 i near  differential  equations 
with  all  admissible  complexities  of  the  field  of  forces.  Such  a  choice  is 
not  necessarily  unique  or  may  not  be  possible  at  all.  The  system  which  we 
shall  recommend  has  merely  been  found  to  be  suitable. 

The  new  system  of  coordinates  will  be  MAR,  an  orthogonal  and  right 
handed  set  of  unit  vectors.  N  is  as  previously  defined,  A  and  R  belong 
to  the  local  horizontal  plane.  A  is  the  unit  vector  of  a  so  that  p  is 
the  angle  between  B  and  A.  Since  the  system  is  right  handed,  ^  is  also 
the  angle  between  T  and  R.  Thus 


(7)  T  =  A  sin  p  +  R  cos  p 
B  =  A  cos  p  -  R  sin  p 
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The  recommended  system  of  variables  is 


V  = 


Vl  = 


V  sin  a,  defined  as  the  vertical  component  of  velocity; 

V  cos  a,  defined  as  the  horizontal  component  of  velocity; 


P 

r 

<f> 


and,  as  may  be  easily  verified,  in  this  system  the  vector  acceleration  (6)  becomes 

(8)  r'  =  (v^  -  pii  )N  +  (Vj:^  +  )A  +  (v^p'  +  ~  cos  p  tan  q))R. 

Notice  that  if  the  sum  of  all  forces  F  (Per  unit  mass)  acting  on  the  point 
is  decomposed  into  its  components  Fj^,  F^  and  Fp,and  if  these  components  can  be 
expressed  as  functions  of  the  cited  variables  only, and  not  of  their  derivatives, 
as  it  turns  out  to  be  the  case  in  our  problem,  then  we  have  a  first  order  quasi- 
linear  system  as  desired,  namely 


(^) 


/  _  n 

''v  ”  r 


>  F, 


,  _  ^h^v 
''h  ■  '  r 


+  F. 


P'  =  -  —  cos  p  tan  <P  +  — 

''h 


r  '  =  V 


=  —  sin  P 


Although  (^)  represents  a  system,  generally  integrable  in  closed  form,  its 
integral  cannot  describe  completely  the  movement  of  a  point  in  terms  of  position 
and  velocity,  because  we  have  only  dependent  variables.  Thus,  one  additional 
dependent  variable  must  be  introduced,  for  example  0  in  the  form 


(10) 


O'  =  Jl  COS  P 
r  cos  f 


Land-Air,  Inc 


5 


GaoiMtrical  Relations 

Basidas  tha  variables  involved  in  (9)  and  (10),  we  are  usually  inter¬ 
ested  in  several  other  parameters  and  their  rate  of  variation.  We  shall 
now  define  those  parameters  and  establish  expressions  that  permit  calcula¬ 
ting  them. 

Let  us  denote  by  Q  the  plane  determined  by  the  vectors  r  and  v.  We 
may  call  Q  the  instantaneous  plane  of  movement  and  it  should  not  be  con¬ 
fused  with  the  osculating  plane  to  the  trajectory.  Q  contains  the  vector 
A  and  thus,  Q  is  also  the  plane  NA.  We  know  from  classical  mechanics 
that,  in  a  central  inverse  square  force  field  (Keplerian  field)  Q  is  an 
invariant  and  later,  we  shall  verify  this  fact  analytically.  We  are  inter¬ 
ested  in  the  angles  determining  the  position  of  Q.  since  their  variation  may 
be  considered  a  measure  of  the  departure  from  Keplerian  movement. 

Let  i,  J,  k  be  the  unit  vectors  in  the  x,  y,  7  directions  respectively 
and  let  r]  be  the  angle  between  Q  and  the  x,  y  plane,  positive  if  |p|  < 

Then 


cos  r]  =  K  •  R 


But 


R  =  T  cos  p  -  B  sin  p 


Hence 


cos  T]  =  K  •  T  cos  p  -  K  •  B  sin  p 


Now 


K  •  T  =  cos  9  and  K  •  B  -  0 


Hence 


(11)  cos  r)  =  cos  <p  cos  p. 
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Let  q  be  the  unit  vector  in  the  direction  of  the  intersection  between  Q 
and  the  xy  plane,  such  that 

-  _  K  X  R 
^  ~  sfn  t)  • 

It  follows  that 

(12)  q  sin  q  =  i(sin  9  sin  (p  cos  P  f  cos  9  sin  P) 

+  j(-cos  9  sin  «p  cos  p+  sin  9  sin  P). 


Let  p  be  the  angle  between  i  and  q,  positive  in  the  xy  direction.  Then 
cos  p  =  i  •  q  and  thus 

cos  p  =  1^'  ■  '  (sin  9  sin  (p  cos  p  +  cos  9  sin  P). 

sin  q 

Similarly,  sin  p  *  j  •  q  and  thus 

sin  p  =  1.^ ■  ■  (-  cos  9  sin  (p  cos  P  +  sin  9  sin  p). 

sin  q 

Multiplying  these  expressions  respectively  by  cos  9  and  sin  9  and  adding, 
we  obtain 


sin  q  cos  p  cos  9  +  sin  q  sin  p  sin  9  =  sin  p 


(13)  sin  p  =  sin  q  cos(p  -  9), 

Eliminating  p  between  (11)  and  (13)»  *e  obtain 

(14)  tan  ip  =  -  tan  q  sin(p  -  9). 
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The  minus  sign  in  (14)  Tollows  from  the  condition  that  at  p  =  0,  we  must 
have  =  0  -  j  and  (p  =  r). 

Other  angles  of  interest  are* 

Y,  the  angle  between  v  and  B.  Given  by 

(15)  cos  Y  =  cos  a  cos  p 

the  angle  between  q  and  r.  Given  b/ 

(16)  cos  =  cos  CP  cos(ji  -  0) 


positive  in  the  direction  of  increasing  cp. 


Precession  and  Nutation 

The  rate  of  variation  of  p  is  uniformly  defined  as  precession  through 
the  literature.  Authors  nevertheless  differ  in  their  definition  of  nutation,  but 
we  shall  denote  by  this  name  the  variation  of  f).  In  effect  and  as  will  be  shown, 
the  quantities  p'  and  r\'  are  not  independent. 

Differentiating  (13)  and  substituting  ro'  with  (1),  we  obtain 


p'  =  _ 

sin  B  cos  <0 


r 


tan  fl)  cos  B 


Differentiating  (13)  and  substituting  G'  with  (1),  we  obtain 

p/  =  cog.Jl.cosiP  »  e}  _  ,  sin  T)  sinjfp  -  OJ 

cos  P  cos  P 

+  —  sin  *1  sin(p  -  9) 
r  cos  V 


Substracting  these  expressions  and  noting  that 
-  Ian  ,  cos  p  =  .2) 

COS  cp 


by  virtue  of  (11)  and  (14)f  we  obtain 
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,^_sin_2 - 

'  'sin  p  cos  <p 


cos  n  cos 


cos 


1^,  ,  sinjLsiilpil  ,  0 


Using  the  identities  already  established,  this  expression  simplifies  to 

(17)  tan(p  -  9)  +  p'  sin  »1  cos  ri*  0 

Expression  (17)  can  also  be  obtained  by  differentiating  (14).  An  alternate 
form  of  (17)  is 

(18)  V  =  p' 


with  cos  ’ll  given  by  (1^). 

The  vector  acceleration,  as  given  by  (8),  can  now  be  expressed  as  a 
function  of  the  precession  and/or  nutation  instead  of  p'.  As  may  be  easily 
verified,  the  corresponding  expressions  are 

=  (v;  -  ^)N  +  (v^  +  ^)A  + 


(10) 


+  either 


linlJL  R 

h  sin  (p 


Using  expressions  (11)  and  (13)»  a  corresponding  quasi-l inear  system 
could  be  wr i tten. 


Keplerian  Movement 

This  is  the  simplest  situation  that  can  be  assumed  and  it  has  been  studied 
to  exhaustion  in  the  literature.  Nevertheless,  it  is  the  only  case  where  a  simpl 
analytical  solution  is  possible  and  therefore,  it  serves  at  least  the  purpose 
of  providing  us  with  a  check  case,  i.e.  a  case  where  the  closeness  of  the 
error  bounds  of  our  numerical  integration  procedure  can  be  verified  numerically, 
a  Therefore,  we  shall  study  the  Keplerian  case  with  this  particular  object  in 
mind. 
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With  reference  to  (9)»  we  heve  in  this  esse 


MQ 

—  where  M  is  the  mass  of  the  earth  and  G  is  the  universal 
constant  of  gravitation. 


Besides  “  0* 

The  system  (9)  becomes  therefore 


r  r* 


V  V. 

v'  =  -  ^  ^ 


(20)  I  P'  =  -  ^  cos  P  tan  9 


r'  =  V 


^  sin  p 

where  MQ  *  3.982  x  10*®  dyne  cm*/gr.  Integration  of  (20)  yields  directly 

P')  ''h=;7r 


where  and  cg  are  arbitrary  constants.  The  latter  expression  can  be  written 
in  the  form 

(22)  V*  =  ^  -f  cg 

Integration  of  (20)  also  yields 


cos  w  cos  P  a  -1.  S  cos  t| 
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where  cs  is  another  arbitrary  constant,  thus  showing  that  the  angle  f)  is  constant. 
This  same  result  can  also  ba  obtained  from  (19).  In  effect,  since  Fp  *  0, 
we  have  f)'  -  =  0  showing  that  not  only  (]  but  also  (x  is  constant,  i.e.  the 

Q  plane  is  invariant.  Since  quite  obviously,  v^  =  t|)'r,  we  have  by  virtue  of  (21) 


(23) 


where,  by  virtue  of  (11),  (13)  and  (16),  vj)  can  be  obtained  from 


I  _  tan  P 
cos  tb  ~  I'  '  r 
^  tan  q 

Let  r'  “  ^  ’1’^  “  ’^‘pc  This  transforms  (22)  into 

(25)  ^)*  -  '’^(cacjr*  +  2cJ  MGr  -  1) 


and  the  general  integral  of  (25)  is 
^  ^  T  -  e  cos?tr+"“c4’) 

where  C4  is  another  integration  constant.  (26)  is  the  equation  of  a  conic 
and  for  |e|  <1,  it  is  the  equation  of  an  ellipse  of  eccentricity  e.  The 
values  of  the  parameters  e  and  p  are  given  by 

*  =  >/^  TSlfisT  ' 

and 

where 

a  =  semi-mtajor  axis  of  ellipse  “  “  ^  • 
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It  follows  that  C8  nHjst  be  negative  for  trajectories  designed  to  remain 
in  the  gravitational  field  of  the  earth.  In  equation  (26),  the  semi-major 
axis  will  be  characterized  by  t|>  =  -  C4. 

In  order  to  obtain  a  relationship  with  time,  we  integrate  (22)  directly, 
obtaining 

(29)  ^  =  ri  /»»’•*  *  2M6''  -  ^  “’•esi"  *»«...■,  It  Cn 

^  ^  7  (MS)>  +  f|  j 

where  C5  is  a  constant. 

We  now  have  all  formulae  necessary  to  verify  the  numerical  integration  of 
(20).  Given  the  initial  values  of  v^,  v^,  P,  r,  q>  and  t,  the  constants  of 
integration  can  be  obtained  as  lOllows:  Cj^  is  obtained  from  (21);  cg  from  (22); 
cg  and  thus  also  q  from  the  formula  following  (22).  The  initial  value  of  t]) 
can  be  calculated  from  (24),  and  C4  from  (26).  Once  c^  and  cg  are  known,  Cg 
can  be  obtained  from  (29)  and  the  conventions  that  must  be  observed  for  this 
determination  may  be  as  follows:  the  positive  sign  is  taken  for  t  and,  assuming 
that  the  initial  data  corresponds  to  a  situation  prior  to  impact  (missile  case), 
the  initial  argsin  is  chosen  between  -  j  and  0.  In  effect,  cgr  +  M6  is  always 
negative  during  normal  missile  flight.  This  can  be  seen  by  noting  that,  at  the 
apogee 


r 


IZ 


0 


ca 


fol lowing 


so  that 


r^cg  +  MG  =  ^  <  0 

and  furthermore,  rcg  MG  is  zero  at  r  »  a,  a  point  usually  below  the  surface 
of  the  earth.  At  impact,  argsin  will  be  chosen  between  -  n  and  - 
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Regarding  the  determination  of  ^  from  (26),  notice  that  prior  to  apogee, 

'll  <  -  C4  and  after  apogee,  'll  >  -  C4, 

Several  formulae  should  be  recalled  for  the  Keplerian  case.  Given  the  values 
of  r  and  v  at  any  point,  the  semi-major  axis  a  can  be  obtained  from 


(30) 


a  = 


2  - 


rv“ 

w 


|f  in  addition,  the  angle  a  is  known  at  that  point,  the  eccentricity  e  can 
be  obtained  from 


X  -2  _ /.  rv*\a  .  r*v*  sin* 

(31)  .  -^1  -fg-J'* — isr- 


Let  the  index  i  refer  to  initial  conditions  and  f  to  final  conditions. 
Then,  the  angle  A'b  =  'lij.  -  'bj  covered  during  flight,  can  be  obtained  from 


If  the  apogee  occurs  at  'll  =  't  ,  then  'll  -  'll.  can  be  obtained  from 

d  3  1 


be  that  value  of  a.  which  maximizes  A'|i  for  a  civen  value  of  v  .  and 

I  *  VI 

f 

2v» 


Let  a. 
im 

known  values  of  r,  and  r^.  This  would  correspond  to  the  minimum  energy  trajectory. 

a.  can  be  obtained  from 
im 


(34) 


sin*  o. 


im 


-IB.  /. 

■  —  V  ■  wj 
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If  (implying  “  'S'j  =  ©duces  to 


(35) 


r  V? 

1  I 


2  - 


MG 


cos' 


im 


2  - 


r  v? 

1  I 

135” 


a_ 
r  . 


Let  e  =  e  if  a.  =  a.  ,  Then,  equation  (31)  can  be  transformed  into 

fn  1  ifT) 


2r  .r. 
.2  =  _ if 


(3M  -  ’ 

r  .V* 

Define  y  =  (35)  becomes 


(37) 


sin*  0.  *  1  ~  ^ 

im  2  -  y 


cos*  o . 


im  ”  2  -  y 
and  (31)  can  be  transformed  into 

(38)  e*  =  1  -  y(1  -  y)  cos*  a. 

so  that  (36)  becomes 


(39) 


e*  =  1  - 
m 


r^v 


f*f 


Define  z  =  HSE”  ’  expression  (32)  can  be  transformed  into 


U 
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(40)  sin  A\k  =  ^  [  y  sin  a.  cos  a.  (1-7  cos^  Or) 
^  e  '  j  i  I 


f 

where  7  =  2-  “(2-y),  and 


+  7  sin  Oj.  cos  (1  -  y  cos*  aj)] 


=  X  -I 


cos*  Of  =  *  7*  cos*  a • 

T  7  I  ^  1 


In  particjlar,  if  r.  =  r-,  then  O;  -  and  (40)  becomes 

Q  1  I  IT 

(41)  sin  A\|r  =  ^  y  sin  a-  cos  a-  (1  -  y  cos*  a  ) 

dll  1 

and,  if  in  (41)  we  set  aj  =  we  obtain 

(42)  sin  A\)^  = 


We  consider  formula  (42)  quite  interesting  for  estimating  the  capa¬ 
bilities  of  a  given  missile,  since  it  gives  the  maximum  geodetic  coverage 
on  a  spherical  globe.  Notice  the  following  table  of  values 


y  ^ 

0  0 

—i _  2 

1  ‘ 

1  TT 


It  should  be  remembered  that  these  formulae  refer  to  an  inertial 
system,  i.e.  do  not  take  into  consideration  the  rotation  of  the  earth. 
This  means  that  the  true  geodetic  coverage  will  be  larger  than  A\^  if  the 
movement  of  the  missile  has  a  positive  West  component  and  smaller  than 
A\)r  if  there  is  a  positive  East  component.  Thus,  in  case  of  war,  it  is 
advantagous  to  have  our  enemy  on  the  West  side.  The  correction  due  to 
the  earth's  rotation,  is  not  difficult  to  calculate.  For  Keplerian 
movement  as  wall  as  for  any  other  case,  we  determine  impact  latitude  and 
longitude,  as  well  as  flight  time  t  in  the  inertial  system  and  correct 
longitude  only  by  adding  cot,  where  O)  is  the  rotation  of  the  earth,  posi¬ 
tive  westward  and  negative  eastward. 
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Similar  considerations  apply  to  the  determination  of  the  vector  v.. 

As  measured  by  an  observer  rotating  with  the  earth,  a  component  (dR  pointed 
eastward  must  be  added,  which  changes  the  value  and  the  direction  of  v.. 

A  further  refinement  which  must  not  be  omitted,  is  that  the  local  plumb 
line  is  not  usually  perpendicular  to  the  local  horizontal,  as  defined  in 
this  paper. 

This  concludes  our  discussion  and  recommended  formulae  for  the 
Kepler ian  case.  Although  far  from  complete,  this  paragraph  remains 
within  its  objective  of  providing  a  checking  case  and  approximate  initial 
estimates. 


ie  designate  by  this  name  a  particular  form  of  equation  (9)  or  their 
equivalents  in  terms  of  precession  and/or  nutation  for  some  suitable 
expression  of  the  force  F.  To  this  effect,  we  decompose  F  as  follows: 

F^  =  the  force  due  to  the  potential  of  the  earth,  i.e.  gravity; 

F  ~  atmospheric  drag.  This  force  may  be  computed  either  considerir^  or 
neglecting  the  effect  of  earth  rotation,  the  decision  depending  on 
some  numerical  analysis  which  is  pending.  By  effect  of  earth  rota¬ 
tion,  we  mean  that  effect  on  air  drag  only,  i.e.  the  latitude  correc¬ 
tion  indicated  in  the  previous  paragraph  is  not  being  referred  to  and 
is  always  assumed  to  be  necessary  after  the  integration  of  the  equa¬ 
tions  has  been  completed; 

F  =  thrust; 

3 

F  =  extraneous  forces  such  as  the  potential  of  bodies  other  than  the 

4 

earth,  light  pressure,  etc.  In  the  case  of  near-earth  missiles,  this 
force  is  considered  negligible  and  will  therefore  not  be  analyzed  in 
this  paper.  It  should  be  noted  that  these  forces  are  definitely  not 
negligible  in  the  case  of  orbiting  satellites. 

Regarding  F^,  there  are  three  cases  to  be  considered,  namely: 

I.  The  earth  is  assumed  spherical  and  the  mass  distribution  is  symmet¬ 
rical  with  respect  to  its  center  0.  In  this  case,  F  is  col  inear  with  N; 
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II.  The  earth  is  e  body  of  revolution  tnasswise  around  z,  which  is  more  re¬ 
strictive  than  the  same  assumption  regarding  shape  only^  and  thus  may 
be  assumed  to  be  an  oblate  spheriod.  In  this  case,  has  components  in 
the  N  and  T  directions  only; 

III.  The  earth  is  not  a  body  of  revolution,  at  least  masswise.  In  this  case 

has  components  in  the  N,  T  and  B  directions,  in  probable  order  of  de¬ 
creasing  magnitude. 

Case  I  is  actually  a  particular  case  of  II  and,  if  no  other  effects  are  bai 
considered,  corresponds  to  Kepler ian  movement.  Therefore,  it  is  not  considered 
sufficiently  important  to  be  studied  separately.  Although  the  refinement  which 
Case  III  represents  over  case  II  cannot,  a  priori,  be  rejected  as  negligible, 
this  case  has  not  been  well  studied  as  to  date.  Blitzer  has  recently  published 
a  series  of  formulae  for  the  earth  potantial  assuming  axial  assymmetry  but  as 
to  date,  their  experimental  verification  is  not  complete.  Thus,  we  shall 
confine  our  study  to  case  II. 

For  the  potential  of  the  earth,  we  take 

U  *  -  ^  M  +  ^  A  -  sin*  <p1  ^  (5  sin*  q>  -  3  sin  q>) 

(43)  ^  ' 

where 


Aa  =  a*  X  1.620?^  x  I0”* 

A3  =  a*  X  (2.20  ±  .08)  x  lO"** 

A4  =  a*  X  2.1  X  10”® 

s^  -  equatorial  radius  of  the  earth  s  6380  km  and  MQ  has  been  given 
following  formula  (20). 

Formula  (43)  represents  the  first  four  terms  of  the  familiar  expanaion 
of  the  earth  potential  in  spherical  harmonics.  This  expansion,  although  con¬ 
vergent,  does  not  converge  to  the  true  value  because  of  cons ider at iom previously 
indicated  under  case  III.  Four  terms  is  the  maximum  which  may  be  taken  under 
any  circumstance  and  it  is  questionable  whether  the  third  and  fourth  terms  are 
reasonable  at  all. 
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The  sign  of  U  has  been  taken  such  that 


(44)  Fi  =  -  grad  U 

so  that,  in  the  NBT  system,  the  components  of  Fj  will  be 

-  "  7  0»  •f'  order. 

Thus  in  the  NAR  system,  the  components  of  F^  will  be 

_  ^  _  sin  p  ^  _  cos  p  flu 

l!r’  ”  r  7lf*  ~  r  ^ 

in  that  order.  Hence,  in  equation  (0),  we  have 

V  -  au 

/ic\  c  Sin  p  flu 

(45)  .  Fi^ - ^ 

-  _  cos  p  au 
'■  IR  “  “  r  ^ 

Notice  that  if  only  F^  is  being  considered,  then  by  virtue  of  (19)  the 
precession  p  can  be  obtained  from 


and,  if  we  consider  only  the  first  two  terms  of  (43) f  then 

^  =  M8  ^  sio  a, 

so  that 

(47)  p'  =  2MGA, 


cos  q  sin*  0 
vTrTsTnTq 
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We  shall  study  Fg  taking  into  consideration  the  rotation  of  the  earth  as 
explained  at  the  beginning  of  this  paragraph.  The  formula,  generally 
accepted  in  the  literature  and  in  the  Transit  project,  is 

(48)  »(h)v,;, 

where  is  a  dimensionless  coefficient; 

A  -  Vehicle's  cross  section  in  the  direction  v  ; 
s  a’ 

<J(h)  =  air  density  as  a  function  of  altitude  h; 

Vg  -  velocity  of  the  air,  at  and  with  respect  to  P,  assuming  that 

the  atmosphere  rotates  rigidly  with  the  earth. 

By  virtue  of  the  previous  definition,  the  vector  v  will  be  given  by 

d 

(4*?)  V  =  (jr  cos  q)  •  B  -  V 

O 

where  u  is  the  angular  velocity  of  rotation  of  the  earth,  i.e.  assuming  one 
revolution  in  approximately  24  hours,  we  have 

u  s  ,0041666667  degrees/sec. 

Now,  since 

V  =  v^  N  +  v^  cos  pB+v^sinpT 


we  have 


v^  =  -  v^  N  +  (wr  cos  q*  -  v^  cos  p)  B  -  v^^  sin  p  T 

(50) 

=  -  v^  N  +  (ur  cos  V  cos  p  -  Vj^)  A  -  <*ir  cos  ?  R 

and 


V*  =  V*  (v^  -  ur  cos  9  cos  p)*’  +  (ur  cos  9)* 
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Using  (48)  and  (50),  we  obtain 


(51) 


''•N  “  ■  5  ®D  ’.''v 

^•A  “  J  ^8  ''a  °®®  <P  cos  P  -  v^) 

Fap  »  -  5  Cd  Ag  <J(h)  v.  cos  q) 


If  wa  want  to  neglect  the  effect  of  earth  rotation,  we  set  v  -  v  and 

d 

u  *  0, 

Notice  that  if  and  Ft  are  being  considered,  formula  (46)  for  the  pre¬ 
cession  is  not  longer  valid.  Instead,  we  must  use 

*  >;;•  S'jy'n  ^  “  5  ^  ''a  ’) 


In  first  approximation,  A  may  be  considered  a  constant.  Otherwise,  it 

s 

can  be  expressed  as  a  function  of  the  angle  between  v  and  v  and  in  this  form, 
will  not  upset  the  quasi-l inear ity  of  the  equations.  o(h)  must  be  known 
eispiricelly«  preferable  in  analytical  form. 

Regarding  the  thrust  Fs,  it  will  be  considered  as  followst 
The  thrust  F3  will  be  sssuswd  constant  and  the  mass  m  of  the  vehicle 
will  be  assumed  to  vary  linearly  according  to 

m  «  m^  -  (m^  -  m^)  r~^T“ 

*  o 

where  m  *  initial  moss  of  vehicle,  at  time  t  »  t  ; 

0  o 

m^  *  moss  of  vehicle  at  burnout,  t  *  t^. 

Hence 


4 

r 


(53)  m  *  oj  -  agt 

where  a^  and  at  are  two  positive  constants  given  by 
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If  in  addition,  F3  is  supposed  colinesr  with  v,  then 


Using  the  same  technique,  additional  forces  can  easily  be  introduced.  For 
reasons  of  briefness,  we  do  not  carry  our  analysis  beyond  this  point.  Combining 
(9),  (45),  (5I)  and  (55),  we  obtain  rather  general  aquations  of  movement,  namely 


'i  *  r  -  ^  -  5  <=D  ^  “ 


I  -  ^v^h  sin  P  9U  .  1 
''h  "7 - 


r  S  *  I  ^  ''a  -P  cos  p  -  v^)  +  cos 


Man  f  -  ^  J  Cp  Ag  o(h)  v^  ur  cos  f 


r'  s  V 


<  «  n  •  A 
9  *  ~  sin  P 


where  sin  a  *  cos  a  —  and  m  is  given  by  (53)  and  (54)»  Since  F8j^  *  0, 
expression  (52)  of  the  precession  continues  to  be  valid.  Since  in  general,  the 
angle  0  must  be  known,  (10)  must  be  integrated  together  with  (56).  In  the  case 
of  a  satellite,  the  nodal  regression  can  be  found  very  accurately  by  integrating 
(52)  over  a  period  of  revolution, 
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W*  designate  by  Nodes  the  intersection  between  the  orbit  of  a  satel¬ 
lite  and  the  equatorial  plane.  In  each  revolution  of  the  satellite,  the 
nodes  move  on  the  equatorial  plane  in  a  direction  contrary  to  the  projec¬ 
tion  of  the  satellites  movement  on  the  equatorial  plane,  by  an  amount  Ap, 
called  nodal  regression. 

In  order  to  obtain  ^  under  consideration  of  the  earth's  oblataness, 
air  drag  and  effect  of  earth  rotation  on  the  latter,  expression  (52)  must 
be  integrated  numerically  over  a  period  of  revolution.  If  air  drag  and 
associated  effects  may  be  considered  negligible,  a  second  approximation 
of  ^  can  be  obtained  from  the  numerical  integration  of  (46).  A  third 
approximation  can  be  obtained  by  numerical  integration  of  (47)  and  we 
must  remember  that  as  we  increase  the  degree  of  approximation,  we  are  in 
this  case,  decreasing  accuracy.  A  fourth  approximation  is  obtainable 
analytically  by  assuming  the  hypothesis  conducive  to  (47)  and  in  addition, 
assuming  that  during  a  particular  revolution,  the  departure  frott:  Kepler ian 
movement  is  negligible.  This  calculation  was  performed  for  the  first 
time  by  Leon  Blitzer  and  published  in  the  Journal  of  Applied  Physics,  Vol. 
28,  No.  11,  Nov.  1957.  Blitzer's  results  have  been  widely  used,  in  many 
cases  with  epparent  indiscrimination  as  to  expectable  accuracy. 

In  order  to  emphasize  this  point,  we  present  our  own  proof  of  Blitzer's 
formula.  From  (47)  and  (24)*  we  obtain 


=  2MG  Az 


cos 


in*  'if 


sin 


and  since  v^ 


r 


^  =  2MG  Aa  cos 


Let  be  the  time  between  consecutive  nodal  passages.  As  t  varies 
from  0  to  T^,  varies  between  0  end  approx imetely  2Tr,  actually  to  a 
somewhat  smaller  velue.  Hence 


22 
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All  *1  M-'  dt  *  2MG  Ab  cos  tj 


J°¥- 


dt 


2n 

^  2MG  Ag  cos  tij  |r? 

where  cos  n  is  held  constant  by  virtue  of  the  assumption  that  we  have  Kepler ian 
movement . 

In  order  to  eliminate  i|i't  we  use  expression  (23),  obtaining 


(57) 


Ap  a  2IIIG  Aa  C!  cos  sir^ 


dij) 


c^  can  be  caluclated  by  considering  the  situation  at  apogee.  Here  r  *  a(1  e) 
and  V.  *  V  where  e  is  the  eccentricity  of  the  particular  orbit  and  v  the 
velocity  at  apogee.  (Do  not  confuse  with  velocity  relative  to  aii^.  Hence, 
since  ^  ,  we  have 


c,  a 


V  is  obtained  from  (30)  by  setting  r  *  a(1  +  e).  We  obtain 

3 


s  /ffiUZI 

“"y  9  I  +  e 


so  that 


(58) 


Furthermore,  by  virtue  of  (26)  and  (28) 


r  -  >(1  -  •*) 

“  1  -  e  cos  (t^  ♦  C4) 


Substituting  this  expres  <on  and  (58)  in  (57),  we  obtain 
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2ii 


J  sin*  tJ)  [l  -  e  cos  (\|)  +  c*)]  d\|) 

tm  sin*  ij)  cos  (ili  +  C4)  dtji  *  0  for  all  constant  04  (another  approx i mat  ion)  and 

Jsin*  \1)  di|)  *  n, 
o 


Now 

2it 


^o 
Hence 


(59)  is  Blitzerte  formula. 

One  of  the  objects  of  this  paper  is  to  calculate  the  departure  of  (59) 
from  the  numerical  integrals  of  (52),  (46)  and  (47), 

Advance  of  the  Line  of  Apsides 

The  line  of  apsides  is  defined  as  the  radius  vector  corresponding  to  apogee. 
Let  'll  1  and  'll  •  be  the  values  of  'll  ,  the  value  of  'll  corresponding  to  apogee  - 
in  two  consecutive  orbital  revolutions  of  a  satellite  and  in  that  order. 

The  advance  of  the  line  of  apsides  is  defined  as 

(60)  A'lig  a  'lig8  -  ^ai 

and  this  movement  is  called  retrograde  if  A'li^  is  negative. 

Again,  we  list  approximations  in  order  of  decreasing  accuracy.  A  first 
approximation  of  A'li^  can  be  obtained  from  the  numerical  integration  of  (56). 

Let  qi  and  P  be  the  values  of  qi  and  p  at  r  ~  maximum.  Then  by  virtue  of  (11), 

d  d 

cos  *  cos  ^  cos  P 

d  d  d 

tan  P 

\  •  Tsrr 

d 


so  that  'll  can  be  calculated  from 

d 


sin  'll. 


s^ 

sin  t) 


ani 


d/or 
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by  virtue  of  (24)*  A  similar  value  of  is  obtained  at  the  next  maximum  of  r 

O 

and  substitution  in  (60)  yields  Aiji  . 

d 

A  second  approximation  can  be  obtained  by  neglecting  in  (56),  all  forces 
not  due  to  the  potential  of  the  earth,  A  third  and  important  approximation 
corresponds  to  the  previous  case  where  we  consider  for  the  earth's  potential, 
only  the  first  two  terms  of  its  expansion  in  spherical  harmonics  as  given  by 
(43)*  The  system  (56)  becomes,  in  this  case 


(61) 


-m 


V  V, 
V  h 


*  -  -  MG  3  sin  2<p  sin  p 


P'  * 


h  A 

-  cos  P  tan  9  -  MG  :  sin  2(p  cos  p 

r  v^r* 


*  ■—  sin  B 


By  using  the  angles  ^  and  f)  instead  of  (p  and  p,  the  system  (61)  can  be 
replaced  by  the  equivalent  system 


(62) 


V  V 

a - -  ^  n  sin  2'|i 


57?fr-  sin  2t|  sin  2^^ 
h 


Details  of  the  transformation  are  rather  obvious  and  are  being  omitted  for 
reasons  of  briefness.  For  the  calculation  of  the  line  of  apsides,  the  system 
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(62)  should  be  preferred  over  (61 )«  since  its  numerical  integration  is 
somewhat  faster  and  thus,  saves  machine  time. 

A  fourth  approximation  can  be  obtained  from  the  analytical  inte¬ 
gration  of  (62)f  a  process  to  which  we  shall  refer  again  with  somewhat  more 
detail,  in  a  later  paragraph  of  this  paper.  It  yields  a  widely  known 
formula  namely 

(63)  /Sir,  =  - -  (2  -  i  sin* 

generally  also  attributed  to  Blitzer  and  published  in  the  same  reference 
of  the  previous  paragraph.  The  approximations  introduced  by  (63)  and  not 
involved  in  (62)  are  of  two  kinds:  first  a  and  e  are  respectively  the 
semi^najor  axis  and  eccentricity  at  the  first  apogee  wA^ere  the  inclination 
of  the  Q  plane  is  r]  and  must  be  calculated  using  the  assumption  that  the 
movement  is  Kepler ian  at  that  time;  second,  ^  as  given  by  (63)  repre- 
sents  some  average  value  of  the  advance  of  the  line  of  apside,  obtained 
by  neglecting  the  terms  which  indicate  its  dependence  on^jr^^.  This  point 
will  be  clarified  later. 

In  general,  the  same  comments  apply  to  formula  (63)  and  (59)*  Never¬ 
theless,  (63)  illustrates  a  fact  which  can  be  easily  evidenced  by  physical 
reasoning,  namely  that  the  line  of  apsides  advances  if  the  orbit  is 
sufficiently  close  to  the  equator,  has  retrograde  movement  if  the  orbit  is 
sufficiently  near  polar  and  is  stationary  at  some  intermediate  value  of  r). 
This  value,  according  to  (63)«  is  approximately  63.5°  and  is  not  in  good 
agreement  with  experimental  observation. 

Thi  Iniigrtl  9f  i  ^miMl  HlBPiititaitriiangiw 

The  case  that  we  are  referring  to,  corresponds  to  the  systems  (61 )  or 
(62)  of  the  previous  paragraph.  The  results  are  only  partial,  that  is, 
these  systems  wore  not  completely  integrated  and  only  the  values  of  certain 
orbital  parameters  were  obtained.  Originality  for  these  calculations  is 
usually  attributed  to  L.  E.  Cunningham,  although  to  the  best  of  the  writers 
knowledge,  the  derivations  have  not  bean  published  yet  as  this  is  being 
written.  The  results,  in  somewhat  different  form,  can  be  found  in  the 
book  by  T.  E.  Sterna,  An  Introduction  to  Celestial  Mechanics,  pp.  122-124* 
For  reasons  of  briefness,  wa  must  omit  the  derivations. 
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It  is  being  assumed  that  the  trajectory  can  be  decomposed  into  infin¬ 
itesimal  segments,  each  one  representing  a  piece  of  a  Keplerian  orbit.  This 
assumption,  obviously  does  not  represent  an  approximation,  but  it  requires 
that  we  sharpen  the  definitions  of  our  parameters;  namely,  let  a,  e,  rj  and 
>)r  represent  respectively  the  semi-major  axis,  eccentricity,  inclination 
of  Q  plane  and  apogee  angle  for  some  £gr,t^£yj^  position  along  the  tra¬ 
jectory.  In  addition,  K^,  K^,  K^,  K^,  is  a  sot  of  constants,  de¬ 

pending  on  the  orbital  parameters  at  that  particular  point.  The  deter¬ 
mination  of  these  constants  will  be  clarified  later.  The  following  formulae 
will  represent  a,  e,  r],  p.,  and  the  mean  anomaly  m  as  a  function  of  \]r. 

The  mean  anomaly  is  given  by  m  =  nt  where  t  is  time  since  apogee  passage 
and  n  is  the  angular  frequency  given  by  n 

a(\|r)  =  K  {•(1  ^  sin*  q)  cos  (ilr  -  ijf  ) 

+  ^  e(1  +  |-)  sin*  q  cos  ^ ^  sin*  q  cos(>)^^^  -  \|f) 

(64)  +  ^(1  -  ^  sin*  q)  cos  2(\)f  -  ^  (1  +  ^*)sin*  q  cos  Zjr 

^  (1  -  I  sin*  q)  cos  3(\)f  -  +  |e  (1  +  J-)  sin*  q  cos(»r  - 

+  |e*  sin*  q  cos  2(^)^  -  +  ^  sin*  q  cos(5^  - 

•(i|r)  ==  K  ^  {(1  f")  (1  -  4  sin*  q)  cos  ()|f  -  ) 

•  a**  (1  -e*)*  ‘■4  2  a, 

■♦■^(1  sin*  q  cos(\|a  +\|r3^)  ''■fgsin*  q  cosOlf^^  -  \lr) 

(65)  ^(1  -  ^  sin*  q)  cos  2(f  “  4  •  H  ^ 

*  ^(1  -  ^  sin*  q)cos  3 (>1^  -  1^3^)  ^(7  +  ^*)  sin*  q  co8(3i|r  - 

+  ^  sin*  q  cos  2{2^  -  \)f3^)  ^  sin*  q  cos(5\)r  - 
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(66)  nC'l')  *  te  co8(\|)  +  \J)^i)  +  cos  2t|i  +  |  cos (3^  -  'I'gi)] 

(67)  li.('li)  -  ■*“  sm(\l>  -  ^t^^)  -  e  sin(>|)  + 

-  sin  2'y  -  |  sinOrb  -  ''bgi)] 

'•’aW  "*  +  iia2~(^-~e^i  {(4-5  sin*  n)  ('b  -  (2  -  3  sin*  r\) 

+  e*  (j  -  ^  sin*  n)  +  ^  sin*  r[  cos  2\b^j]  sin(^  -  ^b^i) 

-  ^[e*  +  j(1  -  5  «*)  sin*  fj]  sin(>|)  +  ^^j)  +  (1  -  ^  sin*  fj)  sin  2('|>-’^) 

(68) 

-  (1  -  I  sin*  ti)  sin  2^1)  +  ^(1  -  ^  sin*  n)  sin  3(tl>  -  ^^j) 

-  ^[«*  "  5(7  +  «*)  sin*  n]  sin  (3'l>  -  >>^1) 

+  ~  sin*  Tj  sin  2(2'|>  -  ibgi)  +  |  sin*  f)  sin(5'l>  -  3’bgi)} 

m(^)  *  +  n^t  +  eaa(1^-  ea')8/»  H(1  ’  |-)  (^  "I  sin*  n) 

+  ^  sin*  T)  cos  2'|i^i]  sin(^  -  'b^j)  +  ^(1  +  ^  e*)sin*  »i  sin(’b  -  ’b^i) 

(6^?)  -  5(1  -  I  sin*  n)  sin  2(^  -  tb^j)  “  %  (1  “  |  sin*  n)  sin3(tb  -  tji^i) 

-  1^(7  -  J^)  sin*  n  sin(3’b  -  'b^i)  -  |  o  sin*  i\  sin  2(2tb  -  ^^^) 

-  ^  sin*  i\  sin(5^  -  3’b^i)} 
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where,  in  (69),  t  is  the  time  since  apogee  passage  at  m  -  m(T|>)  and  n  is  the 
mean  anomalistic  motion,  namely 


(70)  n^  «  n[l  +  (1-3  sin*  p)  ] 

n  being  the  angular  Trequency  at  the  reference  point,  namely  *  9  r 
corresponding  also  to  the  reference  point.  ^ 

The  constants  K  must  be  determined  so  that,  when  i|i  takes  the  value  corre¬ 
sponding  to  the  reference  point,  a(^)  ~  a,  e(^)  -  e,  etc. 

We  are  particularly  interested  in  the  variation  of  the  parameters  p  and  >t> 

d 

and,  from  (67)  and  (68),  it  appears  that  this  variation  is  not  periodic.  In 
order  to  facilitate  calculations,  we  set  the  reference  at  the  first  apogee. 
Then,  if  t|)  •  i|>  ^  is  approximately  equal  to  2ii,  we  will  be  at  the  second  apogee. 

d 

Hence 

Aa  cos  n 

(71)  (i(+„)  =  -  STii  -  .>J»  [*11-  (1  sin 

Af 

Kg  ~  -  a*)*  -  5  sin*  tj) 

•  -  <2  *1^')  “i"*  lJ  +  1  -  i(13  +f)  sin*  sin 

Neglecting  the  periodic  terms,  in  order  to  obtain  an  average,  we  have 
approximately 


cos 


2x  Aa  CO 
■  .>  (1  -  6 


and 


Now,  if  wo  neglect  the  periodic  terms  in  (67)  and  (68),  use  as  reference 
point  the  first  apogee  and  set  we  see  readily  that  »  l*(’^al) 

and  Hence 
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f73)  A,  .  - 
and 

(74)  *  a’l^f  71)1  (4  -  5  sin*  n) 

(74)  is  identical  to  (63);  the  negative  sign  in  (73)  indicates  that  the  nodal 
precession  is  retrograde^  so  that  (73)  is  equivalent  to  (59)  which  was  derived 
by  a  direct  procedure.  The  formulae  (71)  and  (72)  should  be  considered 
superior  to  (59)  and  (63)  respectively,  if  the  departures  from  the  averages 
are  desired. 

On  the  Integration  of  4uasi -linear  Systems 
Consider  a  vector  variable 

X(X^,  Xg,  X3, ...x^) 

and  its  derivative 

x'(x(,  xg,  Xa*****^^^ 

where  primes  indicate  differentiation  with  respect  to  an  independent  variable  t. 
We  are  concerned  with  presenting*  without  proof  for  reasons  of  briefness,  a 
numerical  integration  algorithm  involving  error  bounds,  for  a  class  of  equations 
reducible  to  the  form 

(75)  »  f(*i  t) 

where  f  is  a  vector  function. 

It  follows  from  (75)  that 

dx  a  f (x,  t)  dt 


30 


Land-Air*  Inc. 


an  expression  which,  for  the  discrete  case,  will  be  given  the  following  alternate 
interpretations: 

Ax  «  f  (x,  t)  At 


and 


£al  1  f  (x  +  Ax,  t  +  At)  At 

and  the  function  f  will  be  required  to  satisfy  the  quite  general  condition 

(76)  Ax  »  [o  f(x,  t)  +  (1  -  o)  f(>t  +  Ax,  t  +  At)]  At 

where  0  ^  o  ^  1  for  all  points  x  of  the  integration  dosMin.  Then,  the  integration 
can  take  the  form 

(77)  ^  *  j[f  (x,  t)  +  f  (x  +  Ax,  t  +  At)]  At 
and  the  corresponding  error  Ac  will  satisfy 

(78)  AC  <  j|[f  (x,  t)  -  f  (x  AX,  t  +  At)]|At 

Since  the  expressions  (77)  and  (78)  do  not  represent  an  explicit  algorithm, 
a  second  order  error,  assumed  negligible,  will  be  introducted  by  using,  in  the 
right  member  of  (77),  the  value  Ax*  f (x,  t)  At.  As  is  obviously  a  vector  in¬ 
crement  and  expression  (78)  should  be  intepreted  in  the  sense  that  the  components 
of  AS  are  half  of  the  absolute  values  of  the  differences  between  the  components 
of  the  two  functions  involved. 

An  integration  routine,  for  arbitrary  f,  should  have  the  following  structure. 
Given 


and  the  functions 
Land-Air,  Inc 


31 


calculate 


(^if*  ^8f »  •••  ^pf)*  ^f»  ^(^i»  ^a» 


subject  to  either 


r 


or 


where  r  and  K  are  known.  All  data  and  results  should  be  printed  out  and  suitable 
entries  should  be  provided  for  the  functions  f  . 

ni 

Procedure. 

1.  Set  xj,  xj,  this  step  is  executed  for  the  first  time,  the 

values  X.  and  t.  should  be  used;  at  every  further  execution,  use  the  results 
from  the  previous  calculation. 

2.  Calculate 


(^i)o  *  t)  M 

(AX8)q  *  fsCx,  t)  At 


« 

using  the  values  from  step  1.  for  x  and  t 
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3.  Calculate 


t)  +  +  (Ax)^,  t  +  At]^  At 

^8  *  5(^8 (x*  t)  +  fa[x  +  (Ax)^,  t  +  At]]  At 


t)  +  f^[x  +  C'^)q»  t  +  At]]  At 

4.  Obtain  the  components  of  Ac  as  follows 

ACi  *  jlfi(x,  t)  -  fi{x  +  Ax,  t  +  At)|  At 
ACa  *  jk8(><»  -  f8(x  +  Ax,  t  +  At)l  At 


^'n  “  "  ^n^^  +  aJ,  t  +  At)l  At 

5.  Calculate  the  cumulative  error  e(e^,  ^8»***‘p)  adding  the  results  of 
step  4*  into  the  corresponding  locations. 

6.  Calculate  new  values  of  the  variables,  in  the  form 


xj  +  Axj 


xa  +  Axa 


t  +  At 

and  return  to  1.  until  x^^.  or  r^  takes  the  specified  value.  If  this  occurs, 
print  out  results. 
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CMrfiUflrtfti  tf  Mifafiiani 

All  our  formulas  refer  to  a  geocentric  inertial  systam  as  described 
in  the  first  paragraph.  The  six  scalars  necessary  to  describe  position 
and  velocity  of  a  pointy  have  been  substituted  in  most  cases  by  the  set  of 
variables 

(79)  r,  v^,  v^,  0.  9,  p 

and  if  other  variables  are  being  used,  the  necessary  transformations  are 
indicated.  Nevertheless,  data  given  and  requested  by  military  organiza¬ 
tions  are  usually  referred  to  a  different  coordinate  system.  The  latter 
will  be  described  in  this  paragraph  and  in  addition,  we  shall  indicate  the 
formulae  to  transform  from  this  military  system  to  our  inertial  system  and 
vice-versa. 

The  military  system  is  right  handed,  orthogonal,  usually  denoted  by  XYZ, 
with  the  origin  fixed  at  some  point  on  the  surface  of  the  earth  and  rotating 
rigidly  with  it.  Let  Pi(^)  be  this  origin,  where  R  is  the  corresponding  radius 
vector  from  the  origin  of  our  geocentric  system.  Then  Z  is  pointing  in  the  ^ 
direction,  i.e.  is  the  local  vertical  by  definition  and  XY  define  the  local 
horizontal  plane,  perpendicular  to  f^.  Y  points  south,  i.e.  belongs  to  the 
Zz  plane  and  X  points  west. 

Let  Rb  be  the  projection  of  R  on  the  xy  plane,  (p  the  angle  R®R  in  that 
direction  and  9^  the  value  of  9  corresponding  to  p  >  R®.  As  before  let  t 
denote  time  and  assume  that  at  t  ^  0,  R®  coincides  with  the  positive  direction 
of  the  X  axis.  Thus 


(80)  0^  »  ut  (u  *  .0041666667  deg/sec) 

9o 

We  may  assume  in  addition  that  at  time  t  ^  the  x  axis  has  longitude  zero. 
Then  9^  is  its  longitude  east  of  P^,  i.e.  -  9^  is  its  longitude  west. 

Let  Xj,  9i,  Zi  denote  the  unit  vectors  corresponding  to  XYZ  and  as  before, 
denote  by  i,  j,  K  the  unit  vectors  in  the  xyz  directions.  As  usual,  primes  will 
indicate  derivation  with  respect  to  time.  Then 


(81) 


a 


7. 

It 
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i  sin  9^  -  J  cos  9 

i  cos  9^  sin  f p  ♦  J  8p  ^o 

i  cos  9  cos  ♦  j  sin  9^  cos  •  K  sin  a 
0  0  0  0  o 
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and 

Xi  -  u(i  cos  0Q  +  j  sin  9^) 

7i  *  u  sin  i  sin  +  j  cos  0^) 

Zi  *  w  cos  ®  (-  i  sin  0  +  j  cos  9  ) 

*  Q  '  O'*  0  ^ 

With  reference  to.  our  inertial  system^  if  the  point  Pp  the  movement  of 
which  is  being  studied,  is  defined  in  the  XYZ  system  by  the  vector  r^  and  if 
R  is  the  magnitude  of  R,  we  have 


(83) 


r  a  RZi  +  r 

0 

*  R(i  cos  9  cos  9  +  j  sin  0  cos  e  +  K  sin  ®  )  +  XXi  +  YYi  +  ZZi 
00  00  o 

s  i  (R  cos  9  cos  e  -t*  X  sin  9  -f  Y  cos  9  sin  ®  -f  Z  cos  9  cos  e  ) 

'  0  ^0  0  0  O  ^0^ 

+  j  (R  sin  0  cos  ffi  -  X  cos  0,  +  Ysin0  sine  ■fZsin9  cos  v  ) 

0  O  O  0  '0  O  *0^ 

+  K(R  sin  9^  -  y  cos  +  Z  sin  9^) 

,  *  xi  yj  4*  iK 


Since  R  must  be  known  and  9^  is  given  by  (80),  (83)  permits  the  calculation 
of  the  inertial  coordinates  xyz  as  a  function  of  XYZ  and  the  corresponding  time. 
Furthermore 


r '  *  RZl  +  r' 
o 

«  RZl  4-  X'Xi  4-  Y'Yi  4-  Z'Zi  4-  XXj  4-  y7J  4-  ZZj 

*  R  u  cos  9  (-  i  sin  6  4*  j  cos  9^)  4*  X^  (i  sin  9^  -  j  cos  9^) 
O  O'*  o  o  o 

4"  Y*  (i  cos  9^  sin  •  4-  j  sin  9  sin  -  K  cos  9) 

0  0  0  0  o 

4'  Z'  (i  cos  9  cos  9  4>  j  sin  9  cos  e  4>  K  sin  •  ) 

0  0  0  0  o 

4*  Xu(i  cos  9^  4-  j  sin  9^)  4*  Yu  sin  9^(-  i  sin  9^  ♦  j  cos  9^) 

4-  Zw  cos  ®  (-  i  sin  9  +  j  cos  9  ) 

O'  O  ^  0 
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(84) 


r'  «  i  (-  Ru  sin  9  cos  sin  0^  +  cos  9  sin  cp  +  cos  0  cos  «_ 

0  0  0  0  0  0 

+  Xu  cos  9  -Yusin0  sincp  -Zusin9  cos  ©  ) 

0  O  ^O  0  ^0  * 

+  j  (Ru  cos  9^  cos  ©  -  X*  cos  9^  +  Y^sin9  sin®  +Z'sin0  cos  o 

^  OO  O  0  0  0^0 

+  Xu  sin  9  +  Yu  cos  0  sin  o  +  Zu  cos  9  cos  ©  ) 

O  0  0  0  0' 

+  K(-  Y'  cos  <p^  +  Z'  sin  9^) 
s  X' i  +  y' j  +  z'K  »  V 


V 

r 


So  that  (84)  gives  the  inertial  components  of  the  velocity.  The  angles  9  and  9 
can  be  obtained  from 


tar,  9  =  J 


and 


tan  ©  «  ■  ■ 

-/x*  +  y* 


where 


X*  +  y*  *  X*  +  R*  cos*  9^  +(Y  sin  9^  +  Z  cos  9^) 


(Y  sin  9^  4*  Z  cos  9^  +  2R  cos  9^) 

V  follows  from 

V 

V  *  N  •  V 
V 

*  (i  cos  9  cos  9  4>  j  sin  9  cos  9  +  K  sin  9)  •  v 

and  v^  can  be  computed  f rom  *  v*  -  v* , 

The  angle  P  can  be  obtained  from 

(85)  cos  p  *  IjSSS-S  (y'x  -  yx^)  cos*  9 

*  \ 

Notice  that  the  angles  9  and  9^  correspond  to  geocentric  latitude.  The 

corresponding  geodetic  latitude  is  •and  •  and  given  approximately  by 

0 


\ 

( 

t 
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(86) 


where  e*  «  ,0067226701  *  oblateness  of  the  earth.  The  corresponding  radius  of 
the  earth  is  approximately 


(87) 


ae 

-  e*  s!n*  ® 


cos  ^ 
cos  q> 


where  a  *  equatorial  radius  of  the  earth  *  2,093  x  10’  ft. 

In  the  determination  of  the  angles  9,  (p  and  p,  the  following  considerations 
apply.  For  9,  we  have 


tan  9  «  J 

1C  3ic 

as  noted.  If  tan  9  is  positive,  then  0  <  9  <  j  if  x,y  >  0,and  n  <  9  <  ^  if 
x^y  <  0,  If  tan  9  is  negative,  then  j  <  9  <  n  if  X  <  0,and  <  0  <  2*  if  x  >  0, 
In  other  words,  a  standard  argtan  routine  can  be  used. 

Regarding  <p,  we  have 


tan  9  *  end  I9I  ^  5 

ypHTyf  < 

Thus,  9  has  the  same  sign  as  z. 

Regarding  P  as  given  by  (85),  we  must  compute  the  sign  of 


(88)  Vj  *  -  x^  cos  9  sin  9  -  y^  sin  9  sin  9  +  z^  cos  9 

where  v^  is  the  projection  of  v  on  f.  If  Vj  is  positive,  then  0  <  p  <  j  if  cos  p 
is  positive,  and  j  <  p  <  x  if  cos  P  is  negative.  If  v^  is  negative,  then  «  <  P  < 
if  cos  P  is  negative  and  ^  P  ^  2x  if  cos  p  is  positive.  Schematically: 


Signum  Vj 


Signum  cos  P 


Bounds  of  P 


+ 


+ 


0  -* 

X  . 


X 

5 


2 


2x 
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The  preceding  theory  permits  thus  the  transformation  from  the  system 


V  h,  X.  Y,  Z.  X',  Y',  V 
(h  *  altitude  of  above  sea  level)  into  the  system 
(89)  r,  v^,  v^,  0,  <p,  p 
The  inverse  problem  would  be  the  determination  of 
X,  Y,  Z,  X',  Y',  !• 

given  (89)  and  t,  h.  From  these  three  scalars,  we  can  determine 

R  c:  h  +  — Js—— — 

-  e*  sin*  ^  *^o 

where  tan  =  (1  -  e*)  tan 

and  9  s  ut 
o 

Next,  we  compute 

X  •=*  r  cos  9  cos  9 
y  =  r  cos  9  sin  9 
z  ~  r  sin  9* 


The  substitution  of  these  values  in  (83)  permits  the  determination  of  X,  Y,  Z* 
Next,  we  comoute 


(90) 


x'  s  V  cos 

V 

y'  *  V  sin 

z'  »  V  sin 

V 

% 


9  cos  9  -  Vj^(sin 
9  cos  9  +  Vj^(cos 
9  +  V.  cos  9  sin 


0  cos  P  +  cos 
0  cos  p  -  sin 

P 


9  sin  9  sin  P) 
9  sin  9  sin  P) 


and  find  X',Y',Z^  by  substitution  in  (84).  Thus,  the  inverse  problem  can 
also  be  considered  solved. 
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If  the  geodetic  distance,  measured  on  the  surface  of  the  earth,  between  two 
points  and  pg  is  desired,  i»e. 

Aadist.  Pi(e^,  Vi),  P8(08,  Va) 

it  can  be  calculated  with 


(91) 


A  *  jCRi  +  Ra)  ••'g  cos  [cos (©a  -  -  uAt)  cos(V8  “  9i)] 


where 


a 


-  e*  sin!* 


cos  ,o 
cos  vi»8  ' 


and 


tan  ®i,8 


tan  •«.« 

1'  ^  e'J* 


9 


and  At  is  the  time  interval  for  moving  from  Pi  to  Pa. 


Atmospheric  Refraction 

In  order  to  measure  the  parameters  involved  in  the  preceding  analysis,  either 
with  Transit  equipment  or  solely  ground  equipment,  it  is  necessary  to  have  some 
criterion  for  the  refraction  effect  of  the  atmosphere.  Due  to  both  periodic 
srtd  secular  variations,  no  great  degree  of  reliability  can  be  expected  from 
any  such  criterion  and  for  all  practical  purposes,  we  recommend  to  use  the 
"MkHiel  Atmosphere”,  derived  by  the  Central  Radio  Propagation  Laboratories  of  the 
National  Bureau  of  Standards  in  1959. 

Let 

n  m  refractive  index 
N  *  refractivity 

related  by 

N  a  (n  -  1)  •  10* 

and  let 

*  refractivity  at  ground  level. 

Then,  it  has  been  established  that 


(«)  M.  -  ♦  4«22J^] 
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where 


T  s  temperature,  degrees  Kelvin  to  ^  degree; 
p  s  atmospheric  pressure  to  nearest  millibar; 

e^  ^  saturation  water  pressure,  to  ^  millibar.  May  be  read  from  tables 
of  temperature  versus  water  pressure; 

RH  »  relative  humidity  to  I5I 


Furthermore,  for  a  height  H  above  tracking  equipment,  it  has  been  established 


that 


(93) 


where 


-  c  H 

N  »  N  e  ® 

0 


“  *"  ir+TFi 

0 


and 


AN 


-  7.32.-'>055’’ 


It  follows,  for  the  average  refrectivity  N  that 
”  N  N 

’  hI  **0*  *  - ®c-n 

^o  e  _  11  ® 


c  H 

e  e 


The  second  term  of  the  last  member  being  usuelly  negligible.  The  re¬ 
fraction  correction  should  not  be  confused  or  superseded  by  the  squint  angle 
correction,  which  considers  the  angle  between  the  RF  axis  and  measuring  (digital) 
axis  of  a  radar  antenna. 
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Excluding  interior  ballistics^  the  preceding  theory  will  permit  the 
accurate  calculation  of  the  trajectory  of  near-earth  missiles  and  satellites. 
At  this  point,  we  cannot  indicate  which  affects  should  be  taken  into  con¬ 
sideration  in  each  particular  case.  Such  a  decision  would  require  the 
following  programs,  the  preparation  of  which  is  pending! 

Two  programs  to  perform  the  coordinate  transformations  indicated  in 
the  text; 

A  subroutine  for  the  integration  of  quasi-linear  first  order  systems. 
Such  a  subroutine  should  not  be  specific  on  the  functions  involved; 

A  subroutine  for  the  analytical  integration  of  the  Keplerian  ease. 

The  main  object  of  this  program  would  be  to  check  the  previous  one. 

Programs  for  the  different  cases  presented  in  the  text,  in  order  to 
evaluate  each  perturbing  affect  separately; 

Programs  for  the  evaluation  of  nodal  regression  and  apaidal  advance, 
for  different  cases. 

An  accurate  connection  with  Transit  equipment  can  be  easily  made, 
using  the  expressions  of  the  previous  paragraph.  Letely,  extensive  studies 
seem  to  be  in  progress,  which  consider  the  departure  of  the  earth's 
potential  from  axial  symmetry.  It  should  be  noticed  that  any  new  expression 
for  the  potential  which  may  result  from  these  studies,  may  easily  be  treated 
using  the  same  techniques  of  this  paper. 
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